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Abstract 

We consider the viscous motion of a thin, axisymmetric column of 
fluid with a free surface. A one-dimensional equation of motion for the 
velocity and the radius is derived from the Navier-Stokes equation. We 
compare with recent experiments on the breakup of a liquid jet and 
on the bifurcation of a drop suspended from an orifice. The equations 
form singularities as the fluid neck is pinching off. The nature of the 
singularities is investigated in detail. 
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1 Introduction 



A problem fundamental to the study of nonlinear partial differential equa- 
tions (PDE's) is the nature of their singularities. Perhaps the most famous 
(and unsolved) problem is the suspected blow-up of the derivatives of the 
velocity field in the three-dimensional Euler equation [0. Shocks, i.e., dis- 
continuities in the velocity, are the type of singularities displayed by the 



one-dimensional inviscid Burgers equation [24 



Still a different type of singularity has to be expected from three-dimensional 
free surface flow, which we will consider here. A similar study of a two- 
dimensional flow has been conducted recently Surface tension will tend 
to make the surface as small as possible by reducing the radius. The clas- 
sical stability analysis of an infinite cylinder of fluid by Rayleigh [22] shows 
that the radius does not decrease uniformly: Due to the constraint of mass 
conservation the fastest growing mode is the one with wavelength A ~ 9ro, 
where r$ is the radius of the cylinder. Consequently, the fluid cylinder will 
decay into drops of roughly that size. 

Once the radius becomes zero locally, i.e. the original column of fluid 
separates, the description in terms of a radius function breaks down. Hence 
the equations must develop a singularity at that point. Although linear sta- 
bility analysis gives a reasonable estimate of the size of the droplets formed, 
it completely fails to predict the shape of the surface once an appreciable 
deformation of the original cylinder is reached |J. For example it does not 
explain the fact that the cylinder does not break up uniformly. Rather, 
regular size drops are, under most circumstances, followed by much smaller 
"satellite drops" . Even higher order perturbation theory || || gives only a 
qualitative prediction of the unequal drop sizes, but is not able to describe 
the shape of the fluid anywhere close to pinch-off. This is not very surprising, 
because the characteristic time of the linear instability is close to the time 
distance from the singularity, where expansions in the radius and the velocity 
are bound to break down. 

Therefore a complete treatment of the nonlinearities is needed. The full 
Navier-Stokes equation with free boundary conditions is extremely compli- 
cated, for both analytical and numerical studies. The only simulations of 
axisymmetric drops we are aware of were restricted to irrotational, inviscid 
flow [|IJ. But even with this restriction simulations close to the singularity 
become extremely costly, since the neck region requires high resolution. A 
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reduction of the problem to one dimension will give huge savings in computer 
time, making closeups of the singularity possible. 

There already exist some one- dimensional equations for axisymmetric free 
surface flow |15|, 0. Lee [15[] only considers the inviscid case. Bogy's equations 
[§] allow for dissipation but are very complicated in structure and do not have 
a clear connection with the original Navier-Stokes equation. We will therefore 
derive a set of one-dimensional equations by expanding the radial variable in 
a Taylor series and keeping only the lowest order terms of the Navier-Stokes 
equation. Several invariances and conservation laws of the Navier-Stokes 
equation are preserved. This will be the subject of the second section, along 
with a linear stability analysis. 

Integrating the equations near the singularity proves to be very difficult, 
since the problem becomes very stiff due to the large range of length scales 
in the problem. We develop a fully implicit centered difference method. This 
scheme is then modified to treat the convection term vv z by an upwinding 
technique which ensures negative definiteness of the numerical dissipation. 
The numerical scheme is detailed in section three. 

There is a fair amount of work applying one-dimensional equations to 



the breakup of jets [0, |2j, liquid bridges [[TSj or hanging drops ||. There 
is also work in this spirit on films lining a cylindrical tube [|ll|]. Yet a de- 
tailed comparison between experiments and one- dimensional models within 
the nonlinear regime is missing. Therefore, we try to compare experimental 
drop profiles with simulations close to the breakup point. This is found in 
section four for two recent experiments. 

The first experiment H examines the decay of a free jet of water, the 
second observes how a hanging drop detaches after it is adiabatically filled 



out of an orifice pi] . We also produce an example with a high viscosity 
fluid. Simulations and experiments agree very well, giving ample support 
to the idea that droplet breakup can be well described by one-dimensional 
equations. 

In the next section we take a closer look at the pinch region. We discuss 
a similarity theory for the nonviscous case and explain its failure. All viscous 
solutions are determined by universal scaling functions close to the pinch 
point. The concluding section briefly discusses the approximation used in 
relation to other types of approximations, its higher order versions, and the 
full Navier-Stokes or Euler equations. Finally, we indicate directions of future 
research. 
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2 The equations of motion 



We start from the Navier-Stokes equation for an axisymmetric column of fluid 
with kinematic viscosity u, density p, and surface tension 7. In cylindrical 
coordinates it reads ]14| 



d t v r + v r d r v r + v z d z v r = —d r p/p + v(c%v r + d 2 z v r + d r v r /r — v r /r 2 ), (1) 



d t v z + v r d r v z + v z d s v z = -d z p/p + v{d 2 T v z + d 2 v z + d T v z /r) - g, (2) 

where v z is the velocity along the axis, v r the velocity in the radial direc- 
tion, and p the pressure. The acceleration of gravity g points in negative 
z-direction. The continuity equation reads 

d r v r + d z v z + v r /r = 0. (3) 

The equations (0) and (|3p hold for < r < h(z,t). The balance of normal 
forces gives 

nan = 7 (i/ J R 1 + l/ J R 2 ), (4) 

where a is the stress tensor, n the outward normal, and Ri and R2 are the 
principal radii of curvature. The tangential force balance is 

nat = 0. (5) 

Explicitly, this gives 

P/p-T^T^\- d r v r + (9zV z )h' 2 -(d r v z + d z v r )h'} = 1 (1/R 1 + l/R 2 ) \ r=h (6) 
1 + n z p 

for the normal forces, and 



—^—[2{d r v r )h' + {d r v z + d z v r )(l -h' 2 )- 2{d z v z )h'} = 0| r= , (7) 
1 + a 

for the tangential forces. The prime refers to differentiation with respect to 
z. Finally, the surface has to move with the velocity field at the boundary: 
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d t h + v z h! = v r \ T=h . (8) 

Since we are going to look at thin columns of fluid relative to their elon- 
gation, we expand in a Taylor series with respect to r. By symmetry we 
get 

v z (z,r) =v + v 2 r 2 H , (9) 

and (0) is satisfied by choosing v r to be 

v r (z,r) = -v' r/2-v' 2 r 3 /A . (10) 

The pressure is expanded in the same way: 

p(z,r) = po +P2r 2 H . (11) 

We now insert (|9l)- (|Tl~l) into (0), and ©-(§) and solve the equations 
to lowest order in r. In the case of (0) this gives 

d t v + v v' = -p' /p + u(Av 2 + v'q) - g. (12) 

Equation (|T|) is identically satisfied to lowest order. 

Remembering that h! is also of order r we get from @ an expression for 
the pressure po in fll2"|): 

p /p + uv' Q = l(l/R 1 + l/R i ) (13) 
Similarly, (0) gives an expression involving v 2 . 

- v' ti + 2v 2 h - v o h/2 - 2v' ti = 0. (14) 
Equations (|13"D and flUD can be used to eliminate p an d v 2 from (jT2|) giving 

d t v = ~v v' - ^(l/Ri + 1/Ra)' + Su(h 2 v' )'/h 2 - g. (15) 

The surface condition (§) says to lowest order 

d t h = -v ti - v' h/2. (16) 
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The formula for the mean curvature + I/-R2) of a body of revolution 

is known from differential geometry ||. Thus, dropping the index on vq and 
denoting the surface tension contribution of the pressure by p, we finally get 



d t v 
P 



vv z -p z /p + 3u(h 2 v z ) z / h 2 - g, 
1 h zz 



(17) 
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[h(l + hl)2 (1 + ^)5 J 



and 



-vh z -v z h/2. (18) 

Here the index 2 refers to differentiation with respect to z. When solving the 
set of equations fll7D , fll8"D for z G [— £] we impose the boundary conditions 



h(±i,t) = h± 



(19) 



and 



v(±£,t)=v±. (20) 
The set of equations (JITD - (fJO) is going to concern us for the rest of 



this paper. We reiterate that the physical velocity field (^), (|10|) described 
by dl7|) , (p~8f) has both radial and longitudinal components with a nontrivial 
r-dependence. The physical pressure flTT| ) also carries contributions from the 
shear stress. This should be born in mind when we refer to v and p in (|T7|), 
([18]) as "velocity" and "pressure". 

There are two important conservation laws for this simplified system. 
First, mass conservation means 



d t V = Trh 2 v\ 



i 1 



(21) 



V = 7T / h 2 dz. 
Second the sum of the kinetic energy 



7T 



E k in = ~P / h 2 v 2 dz, 



(22) 



(23) 
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and the potential energy 



E pot = 27T7 / hyjl + h\ dz + npg / h 2 zdz (24) 
obeys the balance equation 

dt(E kin + E pot ) = 



V - 7T ( f /iV - ^=r + P h2y - ^pvh 2 v z + pgh 2 vz 



(25) 



So, apart from boundary terms the total energy changes with the rate of 
energy dissipation 

V = Snisp J {hv z ) 2 dz. (26) 

Since D is negative definite, it follows that, without external driving (bound- 
ary terms in ( p5|) ) , the total energy can only decrease. Note that the poten- 
tial energy for the full equations is precisely ( p4|) , so that the exact surfaces of 
static equilibrium are also equilibrium surfaces of the model: they are states 
which minimize E po t Ml - Famous examples are the equilibrium shapes of 
pendant drops |19[ . This was the reason for keeping lower order terms in the 
expression for p: in a consistent expansion by orders of r the expression for 
p simply would have been 

p = n / (l/h-h gz ), 

resulting in a different form of the potential energy. We also note that T> is 
not negative definite for the viscous term as cited by Cram . His term vv zz 
may feed energy into the fluid, which we found to prevent the system from 
reaching an equilibrium state. 

Although of limited applicability in practice, it is instructive to repeat the 
stability analysis for a fluid cylinder in the case of our model. Assume a cylin- 
der of radius r receives a sinusoidal perturbation of wavelength A = 2%/k] 
then 



r(z,t) = ro[l + e(t)cos(kz)], 
v(z,t) = e(t)vosin(kz). 
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Assuming e(t) = eexp(ut), (|T7j ) and ([18]) give to lowest order in e 

T 

ujvq = (k/ro — r k 3 ) — 3vvok 2 

and 

uj = —Vok/2, 

respectively. This leaves us with the dispersion relation 

u? = uj 2 ((r k) 2 -(r k) 4 )/2-3uujk 2 J (27) 

The solution of (WK) is 



u> = Uo\J (kr ) 2 (l - (kr ) 2 )/2 + ? ^(A;r )4 - h^)\{kr,f , (28) 
where 

4 = v 2 ph (29) 
is a viscous length scale. Both the limits of zero viscosity, 



uj = u yJ{kr Q ) 2 {l - (fcr ) 2 )/2 (30) 

and high viscosity 



UJ 
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.1 - (kr ) 2 )/6, (31) 
r 0/ oz/ 

coincide with the exact result [Q if an expansion to lowest order in kr$ is 
made. 

Equation ( p8|) shows that there is an instability for long wavelengths, the 
stability boundary being krQ = 1 independent of v. In the case of a random 
disturbance, however, the relevant quantity is the most unstable or fastestest 
growing mode. In the general case this is 

w~=icrT7fT)- (32) 
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For v = 0, the most unstable wavelength is therefore \ max = 8.89 r instead 
of the exact value of 9.01 tq [|J. In the limit of very high viscosity the infinite 
wavelength perturbation becomes the most unstable one. 

3 Numerical procedure 

The numerical approximations were computed using a rather simple finite 
difference scheme. The spatial mesh is highly nonuniform, graded mesh; its 
refinement is based on the behavior of the computed solution. The time- 
integration method is an adaptive fully implicit ^-weighted scheme. 
Let the space mesh be 

z 1 < z 2 < ■ ■ ■ < z N 

and adopt the following notation: 

Azj = Zi+i — Zi, 

= (zi + z i+1 )/2, 
Az l+ , = z i+ s_- Zi _i_. 

The meshes used were always constrained to satisfy 



2 " Az i+1 ~ 

The solution at each time level is defined by two arrays, {hi}^L and 
{vi}^S[ l ; the quantity hi is the value of the approximate radius h at the 
mesh point z t and the value gives the value of the approximate velocity 
v at the point z i+ i. In describing the discrete equations for a particular 
time step it is convenient to let dv{ and dhi denote the changes in vi and hi, 
respectively, that take place over the step. 

Difference analogs of the v-equation, (|T7|), were written corresponding 
to each point z i+ i and the difference analogs of the /i-equation, (0), were 
written for each zi. The time derivative term was approximated by dvij At 
or dhi/ At, respectively. The relation for p was used to define it at each point 
Zi in terms of h at z^\, z iy and z i+1 , using centered differences for the h z 
term and a second difference for h zz . (Near the bottom of a pendant drop 
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this was changed; see the remarks just after (|3^).) This defines p at each 
time level in terms of h at that level. 

In setting up the difference equations that mimic ([T7]) and ( P~8| ) the spatial 
terms (everything except the time derivative terms) are evaluated using a 
weighted average of the current value and the yet-to-be-computed value. 
These "mid-step" values can be written as Vi + 9dvi and hi + Odhi. With 
9 = 0.5 this gives a second order correct in time difference equation, but we 
used 9 slightly larger than 0.5 (typically 9 = 0.55). Using 9 close to one half 
gives a small first order truncation term (say 10% of the first-order backward 
difference equation). Taking 9 > 0.5 gives smoother discrete solutions than 
9 = 0.5. 

The approximation of the vv z term at z i+ i is done as follows: 

vv z = Vi(v i+ i - Vi_i)//\z i+ i + NVT 

where NVT is the numerical viscosity term that "upwinds" this nonlinear 
convective term. The NVT is structured so that it is an energy dissipation 
term of small size; the usual technique of simply skewing the difference equa- 
tion in the direction that the fluid is coming from does not assure such a 
property. The NVT term that we use is a difference analog of 

-^-(h 2 V(z)v z ) z , 

where u(z) = d v Az and $ is a nonnegative parameter. Specifically, 

AziV i+ i + Azi+iVi 

Vi + l = A Ta 

Az { + Az i+1 
h 2 v(z)v z \ i+ i = h 2 i+1 (v i+1 - d v i+ i 

-1 h 2 u(z)v z \ l+l - h 2 v(z)v z \i 



NVT 



((h i+1 + ^)/2) 2 Azi 

The rest of the f-equation formed as central differences. Note associating p 
with the Zi's gives p z at the z i+ i points. The viscosity term in (|T7|) is very 
similar to the NVT term; the i>-term is just the constant v. 

The /i-equation at Zi has two spatial terms. The first, vh z , is approx- 
imated by Vi (h i+ i — hi-i)/(Azi + A^_i). The second, v z h/2, is approxi- 
mated by 

/ Vi - Vj-i 

\Azi + Azi-i 



hi. 
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In solving the nonlinear difference equations we use Newton's method. 
Many of our simulations used only one Newton step per time step, starting 
from an initial guess based on linear extrapolation from the previous two 
time levels. It is quite easy, and reasonably efficient, to control the time 
step in such a way that one step of Newton's method reduces the error to 
very close to rounding error. It is worthwhile pointing out that even if the 
decision is made to only use one step of Newton's method it is useful to code 
it in general, since observing quadratic convergence of the iteration is a good 
check on whether the linearization has been done correctly. 

4 Comparison with experiment 

The first experiment we consider studies the breakup of a liquid jet |J. Water 
is pumped through a nozzle at high speed to form a liquid column virtually 
unaffected by gravity. A periodic perturbation, whose amplitude and fre- 
quency can be controlled, is applied to the jet as it leaves the nozzle, The 
system is allowed to reach a steady state, in which the jet at a sufficiently 
large distance from the nozzle has completely broken up into droplets. Pho- 
tographs of this stationary configuration are taken. 

We try to model the experiment as closely as possible, but since we can 
only simulate up to the point of the first singularity (due to limitations of 
our current program) we cannot reach the stationary state. Instead, we fix 
h + = h- = r <C C and v + — v _ = V, and over a period of 8 wavelengths 
smoothly turn on a small sinusoidal perturbation to t>_. 

Thus the parameters of the simulation are the length of the jet 2£, its ini- 
tial radius r , the fluid parameters 7/p and u, the speed of the jet V, and the 
amplitude V p and frequency f p of the perturbation. We chose r /2£ = 0.004, 
so the size of the drops is very small compared with the jet length and the pre- 
cise value of this ratio is immaterial. V p was adjusted to make breakup times 
conform with experiment. The remaining dimensionless parameters control- 
ling the problem are A/r , £ v /ro, and the Weber number /3 2 = proV 2 /^. 
Here A = V/f p is the wavelength of the perturbation and £ u the viscous 
length dffi) . Typical values for fluid parameters can be found in Table 1. 
The jet experiments were done with water. 

For the jet, the linearized problem of Section 2 now takes place in a semi- 
infinite geometry, where surface perturbations are prohibited to the left of 
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the nozzle opening [|13|,|T6[. However, for large Weber numbers (239 in the 
present experiment) the growth of unstable modes is just the same as the 
temporal growth of Section 2, translated into space via the jet velocity V. 
Also, the parabolic velocity profile of the nozzle opening has relaxed into a 
plug profile in the relevant region of the jet ||, so we are assuming a constant 
profile right from the opening. 

We follow the simulation up to the first singularity. The resulting profile 
is aligned with a picture of the experimental jet, to make the minima in 
front of the drop which is about to detach coincide. In Figure 1 theoretical 
and experimental profiles are compared for A/ro = 14.57. The case with 
the smallest perturbation is shown ||. Allowing for some blur of the pho- 
tographs, the agreement in the shape of the drop about to form is quite nice. 
Note that the breakup is taking place in a very asymmetric fashion (with 
respect to the breakup point): On the right side the profile is quite steep 
forming a very much rounded drop; on the other side a flat neck formed, 
which will eventually coalesce into a smaller satellite drop. 

An even more direct comparison is possible with an experiment investi- 
gating a dripping tap • As long as the drop is small, it will be suspended 
stably from the orifice. By slowly filling in more liquid, the drop goes through 
a series of stable states, until eventually gravity overcomes surface tension 
and the lower half of the drop falls. Subsequently, a thin neck forms and 
the lower part of the drop detaches. The stability of the drop hanging in 
equilibrium has been the subject of much study in itself fl9fl . The length 
scale controlling this problem is the capillary length 

4 = {i/ P g) l/2 . (33) 

For water, £ c and £ u , the viscous length, are separated by almost five orders of 
magnitude, see Table 1. Hence there is a wide range of physical phenomena 
to explore between the onset of the linear instability and the breakup of the 
drop. 

We will not repeat the stability analysis for our one-dimensional equations 
here, but concentrate on the breakup. The only dimensionless parameters in 
the problem are the ratios £ c /r and £ u /ro, where r is the radius of the orifice. 
They were made to coincide with the experimental values, the working fluid 
being water. 

There are some technical problems involved in simulating the moving 
boundary at the lower end of the drop. We avoid having to use a movable 
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grid by mapping the problem on the unit interval, z/£ = x G [0, 1] where £ 
is the length of the drop which is calculated using 



l{t) = f v (l,s)ds. (34) 

Jto 

By definition, v(l,s) is the velocity of the lower boundary. Care must also 
be taken to calculate the pressure at the endpoint where h z becomes infinite. 
For the values x G [0.9, 1] we calculate p by interpolating h(x) with an even 
fourth-order polynomial. Then all the singularities in the mean curvature 
cancel. 

Figure 2 shows a series of profiles taken at constant time intervals of 
0.4(rQp/7) 1 / 2 . In the experiment, this would correspond to 6.6ms. Given the 
very small time scale it would be very costly to let fluid drip as slowly as is 
possible in experiments. To still let initial oscillations die out, the viscosity 
is set to a very high value initially, and is then reduced to the value of water 
well before the first instability. Fluid is injected at the orifice with speed 

1/2 

0.02 7/(pr ). To the profiles at constant time intervals we add a snapshot 
of the drop as the width of the neck becomes 0.01 r . We also superimpose 
an experimental picture of the drop j21 |, taken at the point of breakup. 



The very good agreement with simulations is especially impressive since 
this was not to be expected from a simple one-dimensional approximation. 
In particular in the lower half of the drop the assumption h ^ £ seems to 
fail, but one must remember that this part of the drop is almost static in a 
moving frame of reference. But the static limit of the equations is retained 
exactly in the approximation. Note that although the linear instability of 
the hanging drop is not investigated explicitly, it is also accurately described 
by the model. Namely, it determines the total volume of the drop (upper 
and lower half combined) and influences the point of breakoff. 

Again, the breakoff occurs very asymmetrically, as was already observed 
in the jet decay. The asymmetry therefore does not come from the action of 



gravity. This is also confirmed by the estimate of Peregrine et al. who 
estimate that by the time the neck is formed, straining forces due to surface 
tension outweigh the straining forces due to gravity. 

We conclude this section by reporting on a simulation of a fluid with 
significantly higher viscosity. With the radius of the orifice being 0.06 cm, 
we adjusted £ c /r and £ u /r to match the parameter values for glycerol at 
25°C, as given in Table 1. The viscosity of glycerol is about 1,000 times higher 
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than that of water, leading to a very different type of dynamics. Figure 3 
shows the neck being pulled into a long and thin thread. Its length is 40 
times the radius of the orifice at the point of rupture. Qualitatively, this is 
consistent with linear stability analysis: for high viscosity, the most unstable 
wavelength becomes large, see equation fl32|) . On the other hand, the radius 
r of the thread becomes very small, so fl3~2|) cannot account for its length 
in any quantitative way. The origin is clearly dynamical. The break occurs 
at the upper end of the thread in the simulation presented, but it may also 
happen close to the drop under slightly different conditions. Experiments 
with high viscosity fluids in the same geometry are in in progress p3]| . 



5 Nature of Singularities 



We now look closer at the point where the fluid neck is pinching off. As 
the radius goes to zero, pressure forces are expected to diverge, and the 
small amount of fluid left in the neck region is pressed out of it even faster. 
Therefore, as h min — > 0, where h min is the minimum radius, the velocity 
and higher derivatives of both h and v will probably become infinite at the 
point of rupture. This is the singularity or "blow up" we want to investigate 
further. 



Keller and Miksis |12| present a very interesting scaling theory for the sin- 
gularity in the nonviscous case. There are two important differences between 
our problem and theirs: Their Geometry is two-dimensional rather than 
three-dimensional- axisymmetric, and they study the time after the breakup. 

The idea of their study may be described as follows: Since h becomes very 
small near the singularity and v large, the pinch region is separated in scale 
from the boundaries. Therefore boundary conditions become irrelevant and 
the flow is determined by 7/p alone. If At — t s — t represents the time dis- 
tance from the singularity, the only available length scale is the combination 
( r yAt 2 /p)3. Hence, introducing 

z = (z-z s )( 1 At 2 /p)--i (35) 
where z s is the position of the pinch point, and 

h = h{^At 2 /p)~^, v = u(7/pAt)~s (36) 
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the problem can be written in terms of the similarity variables z, h, and v 
alone. Once the similarity equation is solved, the resulting profile determines 
the evolution of the interface for all times up to the singularity. Note that 
h — ► and v — > oo as At — > if h and v are assumed fixed, consistent with 
the original assumptions. 

We will see, however, that this similarity argument does not carry through 
for the case of our equations, since the inviscid case appears to develop 
singularities even before h m i n — > ! For a consistent formulation up to 
the point of breakup we therefore need to add at least a small amount of 
viscosity. We are confident that the following conjecture, due to Constantin 
[§], is true. It indicates that with viscosity the singularity does not occur 
until h min goes to zero. 

Conjecture 1 Forv > andt G [0, to] such that h(t) > h > the solutions 
°f vHj)' HHb) s ^ a U regular, i.e., h, v, and all their derivatives remain bounded 
in [—1, 1] for t G [0, to] with bounds depending only on v and h . 

To investigate the problem further, we chose a cylinder of radius r = 0.01 
and length 2 as an initial condition. At its ends z = ±1 the radius is kept 
fixed and the velocity is set to zero. Given a slight initial disturbance in 
the velocity field, the cylinder collapses and forms a singularity after about 
20(r^p/7)i The viscous length is l v = 1.4 ■ 10" 5 . 

The profile near the singularity comes out quite asymmetric, as observed 
before. This is also true if one starts from initial data almost symmetric 
around z = 0. Two almost linear pieces of different slope are joined smoothly 
by a round piece with a radius of curvature comparable to the minimum 
radius, cf. Figure 4 . Hence both terms in the pressure in fli~7|) are of the 
same order of magnitude at the minimum, while the first term dominates the 
linear regime away from the minimum. This means the pressure is higher on 
the shallow side of the minimum (right hand side in Figure 4), forcing fluid 
over to the steeper side (left hand side in Figure 4). As seen by comparing 
with the velocities in Figure 4, the minimum is convected with the velocity 
of the fluid, so in the frame of reference of the minimum fluid is expelled 
on either side. At the same time this causes the left hand side to get even 
steeper. 

If there is no mechanism curbing this process, the slope will eventually 
get infinite. All our simulations, conducted with different initial radii and 
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initial disturbances, show that for the inviscid equations exactly this happens 
and h z goes to infinity even before h m i n goes to zero. Note that t v is still 
less than 2 % of the minimum height in Figure 4, yet the inviscid equations 
already would have blown up at the times shown. If v is finite, on the other 
hand, h z turns out to be uniformly bounded by a constant which gets larger 
as v decreases. Hence for finite, but arbitrarily small v blow-up only occurs 
for h min -> 0. 

From the conjecture mentioned earlier we conclude that viscous solutions, 
up to some finite minimum radius ho, can be well approximated by finite 
differences as long as the mesh is fine enough. Intuitively, one expects that the 
mesh size Az should be at least of the order of ho- We checked convergence 
near the singularity by conducting a series of runs with increasingly fine 
resolution. To save on computational effort, only the region around the 
singularity is highly resolved, the grid getting coarser by factors of 2 towards 
the outside. We plotted h m i n and the maximum velocity v max versus the time 
difference from the singularity on logarithmic scales, see Figure 5 . Lengths 
are shown in units of l v and times in units of the viscous time scale 

t„ = z/ 3 p 2 /7 2 (37) 

The plots agree up to the length scale of the coarser grid. We also monitored 
the highest derivatives in the problem, i. e. p z and v zz . The dashed vertical 
line indicates up to which point they seemed well resolved. As can be seen 
in Figure 5, problems with resolution occur when h m i n is of the order of 
Az, indicated by the horizontal line. Since the numerical viscosity NVT 
as introduced in Section 3 is approximately equal to vAz, convergence for 
increasingly fine grids also demonstrates that it does not introduce artificial 
effects. For the finest resolution the numerical viscosity was less than a tenth 
of the physical viscosity in the center of the grid. From all this we feel 
confident that the plot of our best-resolved run in Figure 5 gives a faithful 
description of the original equations up to the point indicated. 

Figure 5 indicates that v max goes to infinity as the singularity h m i n 
is approached. All derivatives of the velocity as well as second or higher 
derivatives of the height are found to blow up even faster. This means their 
asymptotic value increases faster than a negative power of Az as we increase 
the resolution. The maximum value of h z , however, approaches a constant as 
h m in — > 0. This is an important self-consistency property of our equations: 
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The solutions never approach a situation where the surface parametrization 
is bound to break down. 

There are some regions where h min (t) is close to a power-law, but they 
never extend over more than two decades in length scales. In v max there is 
even less an indication of power-law behavior. The decay of h m i n is always 
faster than the ts power-law predicted by fl35|) , (|36|). 

Considering also the profiles h(z) and v(z) directly, we conclude that (135|) . 
(|36D is clearly not valid, even for h min ^> £ v . The reason may be that {h z ) max 
goes to infinity long before h min goes to zero, hence viscosity is important 
even on scales much larger than £ v . The nature of the singularity of our 
inviscid equations probably is specific to the approximation. For example, 
it could be that the full Euler equations, instead of overturning, produce 
localized regions of high vorticity which our approximation cannot describe. 

The system described so far is determined by the four parameters tq, £, 
■y/p, and v . If it has a solution h(z, t) and v(z, t), the system with parameters 
aro, a£, (a 3 /6 2 )7/p, (a 2 /6)z/ will have the scaled solution 



h a b(z,t) = ah(z/a,t/b), 

Vab = T v(z/a,t/b). 
b 

This is equivalent to saying that up to a rescaling of length and time 
the solution is determined by two dimensionless ratios, r /£ and £„/£, say. 
By the argument presented at the beginning of this section, one expects the 
solution near the singularity to be independent of the dimensions of the initial 
cylinder. Hence all solutions, for t ~ t s and z ~ z s , can be written in the 
universal form 

h{z,t) = £ u h s (±(z- z s )/£ u ,(t s -t)/t u ), (38) 
v{z,t) = ±fv s {±{z-z s )/£ v ,{t s -t)/t u ). 

The ± signs in (|38D take care of the fact that the solutions may have 
different parity, with fluid flowing from left to right or vice versa. 

We tested ( |3~8D by conducting simulations with different parameter values 
and calculating h s and v s from them. Namely, we increased r by a factor 
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of 10 and also varied the viscosity. This causes the global behavior of the 
solution to change dramatically, yet on length and time scales comparable 
with l v and t v or smaller Q58D is obeyed beautifully. Figure 6 shows h s and 
v s as calculated from different runs, all at (t s — t)/t v = 1.97. Note that 
the reduced profiles still evolve in time, unlike the similarity solutions of 



6 Discussion 

The key to the success of the present investigation lies in the construction 
of appropriate model equations to study the motion of thin columns of fluid. 
First, our expansion method allows to take viscous body forces as well as vis- 
cous boundary conditions into account. This makes it distinct from methods 
where the average velocity over the cross section is the dynamical variable, 



such as in the equations for shallow water waves p0| . Precisely due to the 
inclusion of boundary conditions, the viscous terms in our equations become 
purely dissipative. 

Second, we take the exact curvature term (|T7D into account. The impor- 
tance of those higher order terms of the expansion for strong variations of h 



was noticed before [fTI]] . Figure 2, for example, beautifully demonstrates how 
the model takes equilibrium shapes into account. Also, regions of high slope 
(h z fa 10) at the top of the drop are very well represented. 

Apart from experimental test, though, we do not see how to give a pri- 
ori estimates of the quality of approximation in the framework of our model 
alone. A possibility is to study the next order in the expansion. Apart from 
an equation of motion for h, we now have two equations for the expansion co- 
efficients of the velocity field, vo(z, t) and vz{z, t). Those equations, although 
readily written down, are considerably more complicated than (17), fll8D and 



require new numerical methods. Therefore, we consider it a study all of its 
own which should be investigated separately. 

Most importantly, our equations remain self-consistent right up to the 
point of rupture h m i n — > 0. This means no other singularity occurs before 
that point. This is supported by our simulations over a wide range of viscosi- 
ties and by preliminary mathematical analysis ||. Specifically, there is no 
overturning of the profile. The full equations of motion certainly would not 
form singularities even in the case of overturning, but there does not seem to 
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be experimental evidence for this to occur before breakup. (This excludes ini- 
tial or boundary conditions with strong transversal velocity gradients, which 
"force" the flow to overturn, but which are not realizable in our equations 
in the first place.) Hence we see no reason to doubt the applicability of our 
model even for small viscosities such as in water. 

However, the inviscid version of our equations clearly is at odds with ex- 
perimental evidence, showing overturning on experimentally accessible time- 
scales. This reflects the singular nature of the limit v — > 0. We hope this 
will shed some light on the nature of this limit in the Navier-Stokes equation 
and its relation to the Euler equation. 

We plan to develop a code with adaptable grid, which moves with the 
position of the minimum and introduces new grid points when needed. This 
code is expected to be much more effective and to allow us to reach consider- 
ably higher resolution. We hope this will allow us to explore the asymptotic 
regime even more carefully. 

Another expected benefit of the new code is to be able to go beyond the 
first singularity by introducing a new grid point at the pinch. The equations 
will then be integrated from there with new boundary conditions. This will 
allow us to investigate a new range of phenomena, like formation of satellite 
drops, recoiling, etc. 

In conclusion, we have developed a one-dimensional equation for an ax- 
isymmetric thread of fluid. Computed profiles coincide nicely with experi- 
ments. The inviscid equations are inconsistent, leading to singularities even 
before the breakup into drops. All solutions with v > are described by the 
universal scaling functions h s , v s sufficiently close to the singularity. 
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TABLE 1 





Water 
20°C 


Glycerol 
20°C 


Glycerol 

25°C 


v [cm 2 /sec] 


1.00 • 10~ 2 


11.8 


7.6 


7/p [cm 3 /sec 2 


fluid-air interface 


72.9 


50.3 


50.0 


ic = (i/pg) 1/2 






0.273 


0.226 


0.226 


l v = pis 2 /7 [err 




1.38- 10 _B 


2.79 


1.15 


t u = u 3 p 2 h 2 [f 


•ec 




1.91 • 10- 10 


0.652 


0.174 
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Table Captions 



Table 1 



This table contains the physical parameters for water at 20°C and glycerol 
at 20°C and 25°C The values are quoted from |H| . 

The first line contains the kinematic viscosity the second the surface 
tension divided by density 7/p. The remaining three lines contain character- 
istic length and time scales, g is the acceleration of gravity. 
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Figure Captions 



Figure 1 



Comparison between a decaying water jet (upper) and our simulation 
(lower). We processed the original image so as produce a white background. 
The nozzle is to the left. The point where the first drop detaches from 
the experimental jet has been aligned with the corresponding point of our 
simulations. The horizontal scale has been adjusted as well. The parameters 
are A/r = 14.57, £ u /r = 6.54 • 10~ 4 , and [3 2 = 239. The fluid parameters 
[|] differ slightly from the ones quoted in Table 1, due to additives. 



Figure 2 



Simulation of a drop of water suspended from a circular orifice of radius 
r = 0.26cm. This makes the parameters l c /r = 0.992 and l u /r = 4.89- 10~ 6 , 
compare Table 1. The time distance between profiles is OAfryp/'y) 1 / 2 , starting 
from a point where the drop is already falling. We also superimpose a profile 



at the snap-off and the corresponding experimental picture [21 1. There is no 



adjustable parameter in the comparison. To enhance contrast, we erased the 
background in the original photograph. 



Figure 3 



Same as Figure 2, but with the fluid being glycerol at 25° C and r = 
0.0625cm. The parameters are now Z c /r = 3.61 and I v /tq = 18.3 . Note the 
long neck, which is the trademark of high viscosity fluids. 
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Figure 4 



A closeup view of simultaneous radius, velocity, and pressure profiles close 
to pinch-off. The time distance At from the singularity are log 10 (At) = —4.0, 
—4.1, and —4.2, in units of j/p and r^. The viscosity is v = 0.0037. The 
pressure is higher on the right hand side of its peak, pushing fluid to the left. 
The minimum of the radius h moves with the fluid underneath it. 

Figure 5 

The minimum radius h m i n , maximum (absolute) velocity v max and the 
maximum (absolute) slope (h z ) max as functions of the time distance from 
singularity At, in units of the viscous scales £ v , t u , and v u = £ u /t u . The axis 
are logarithmic. The dashed vertical line indicates the point where p z is no 
longer fully resolved. This happens when h m i n reaches the grid size. The 2/3 
- slope would be predicted by a nonviscous similarity theory. 

Figure 6 

The reduced profiles h s and v s as calculated from different parameter 
values (r Q ,£, j/p, v) via fl38l) . The solid line represents (0.01,1,1,0.0037), 
the dotted line (0.1, 1, 1,0.0037), and the long-dashed line (0.01,1,1,0.0074). 
The point of touch-down is shifted to zero in each case, the units are the 
viscous scales. The dotted lines had to be flipped over (- signs in (|38D ) to 
correct for the difference in parity. 



26 



Figure 1: 




27 



Figure 2: 



1 » 




H 






Bf// 



28 



Figure 3: 




29 




-0.93 -0.925 -0.92 



position z 



30 




31 




(z-z s )/l v 



32 



